#include "ManagedLasoProvider.h"

/* Table of constant values */

static integer c__0 = 0;
static integer c__1 = 1;


/* ------------------------------------------------------------------ */

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::sippla_(S_fp
   op, S_fp iovect, integer *n, real *xl, real*xr, integer *np, integer *nperm,
   integer *nop, integer *nmval, real *val, integer *nmvec, real *vec, integer
   *nblock, real *h__, real *hv, real *p, real *q, real *bound, real *d__, real
   *deltal, real *deltar, logical *raritz, real *eps, integer *ierr)
{
  /* System generated locals */
  integer val_dim1, val_offset, vec_dim1, vec_offset, h_dim1, h_offset, hv_dim1,
    hv_offset, p_dim1, p_offset, q_dim1, q_offset, i__1, i__2, i__3, i__4;
  real r__1, r__2;

  /* Builtin functions */
  //double sqrt(doublereal);

  /* Local variables */
  static real hmin, hmax, temp;
  //extern real sdot_(integer *, real *, integer *, real *, integer *);
  static real zero[1];
  //extern real snrm2_(integer *, real *, integer *);
  static integer i__, j, k, l, m;
  //extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
  // integer *), saxpy_(integer *, real *, real *, integer *, real *, 
  // integer *);
  static integer jj, kk;
  //extern /* Subroutine */ int slaeig_(integer *, integer *, integer *, 
  // integer *, real *, real *, integer *, real *, real *, real *, 
  // real *, real *, real *, real *, real *), slager_(integer *, 
  // integer *, integer *, real *, real *, real *);



  /* THIS SUBROUTINE POST PROCESSES THE EIGENVECTORS.  BLOCK MATRIX */
  /* VECTOR PRODUCTS ARE USED TO MINIMIZED THE NUMBER OF CALLS TO OP. */

  /* Parameter adjustments */
  q_dim1 =  *n;
  q_offset = 1+q_dim1 * 1;
  q -= q_offset;
  p_dim1 =  *n;
  p_offset = 1+p_dim1 * 1;
  p -= p_offset;
  hv_dim1 =  *np;
  hv_offset = 1+hv_dim1 * 1;
  hv -= hv_offset;
  h_dim1 =  *np;
  h_offset = 1+h_dim1 * 1;
  h__ -= h_offset;
  val_dim1 =  *nmval;
  val_offset = 1+val_dim1 * 1;
  val -= val_offset;
  vec_dim1 =  *nmvec;
  vec_offset = 1+vec_dim1 * 1;
  vec -= vec_offset;
  --bound;
  --d__;

  /* Function Body */
  zero[0] = (float)0.;

  /* IF RARITZ IS .TRUE.  A FINAL RAYLEIGH-RITZ PROCEDURE IS */
  /* APPLIED TO THE EIGENVECTORS. */

  if (!(*raritz))
  {
    goto L180;
  }

  /* ------------------------------------------------------------------ */

  /* THIS CONSTRUCTS H=Q*AQ, WHERE THE COLUMNS OF Q ARE THE */
  /* APPROXIMATE EIGENVECTORS. */

  i__1 =  *nperm **nperm;
  scopy_(&i__1, zero, &c__0, &h__[h_offset], &c__1);
  m =  *nperm %  *nblock;
  if (m == 0)
  {
    goto L40;
  }
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    scopy_(n, &vec[i__ *vec_dim1 + 1], &c__1, &p[i__ *p_dim1 + 1], &c__1);
    /* L10: */
  }
  (*iovect)(n, &m, &p[p_offset], &m, &c__0);
  (*op)(n, &m, &p[p_offset], &q[q_offset]);
  ++(*nop);
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 =  *nperm;
    for (j = i__; j <= i__2; ++j)
    {
      jj = j - i__ + 1;
      h__[jj + i__ * h_dim1] = sdot_(n, &vec[j *vec_dim1 + 1], &c__1, &q[i__
        *q_dim1 + 1], &c__1);
      /* L20: */
    }
    /* L30: */
  }
  if (*nperm <  *nblock)
  {
    goto L90;
  }
  L40: m +=  *nblock;
  i__1 =  *nperm;
  i__2 =  *nblock;
  for (i__ = m; i__2 < 0 ? i__ >= i__1: i__ <= i__1; i__ += i__2)
  {
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ -  *nblock + j;
      scopy_(n, &vec[l *vec_dim1 + 1], &c__1, &p[j *p_dim1 + 1], &c__1);
      /* L50: */
    }
    (*iovect)(n, nblock, &p[p_offset], &i__, &c__0);
    (*op)(n, nblock, &p[p_offset], &q[q_offset]);
    ++(*nop);
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ -  *nblock + j;
      i__4 =  *nperm;
      for (k = l; k <= i__4; ++k)
      {
        kk = k - l + 1;
        h__[kk + l * h_dim1] = sdot_(n, &vec[k *vec_dim1 + 1], &c__1,  &q[j
          *q_dim1 + 1], &c__1);
        /* L60: */
      }
      /* L70: */
    }
    /* L80: */
  }

  /* THIS COMPUTES THE SPECTRAL DECOMPOSITION OF H. */

  L90: hmin = h__[h_dim1 + 1];
  hmax = h__[h_dim1 + 1];
  slager_(nperm, nperm, &c__1, &h__[h_offset], &hmin, &hmax);
  slaeig_(nperm, nperm, &c__1, nperm, &h__[h_offset], &val[val_offset], nperm,
    &hv[hv_offset], &bound[1], &p[p_offset], &d__[1], &q[q_offset], eps, &hmin,
    &hmax);

  /* THIS COMPUTES THE RITZ VECTORS--THE COLUMNS OF */
  /* Y = QS WHERE S IS THE MATRIX OF EIGENVECTORS OF H. */

  i__2 =  *nperm;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    scopy_(n, zero, &c__0, &vec[i__ *vec_dim1 + 1], &c__1);
    /* L110: */
  }
  m =  *nperm %  *nblock;
  if (m == 0)
  {
    goto L140;
  }
  (*iovect)(n, &m, &p[p_offset], &m, &c__1);
  i__2 = m;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    i__1 =  *nperm;
    for (j = 1; j <= i__1; ++j)
    {
      saxpy_(n, &hv[i__ + j * hv_dim1], &p[i__ *p_dim1 + 1], &c__1, &vec[j
        *vec_dim1 + 1], &c__1);
      /* L120: */
    }
    /* L130: */
  }
  if (*nperm <  *nblock)
  {
    goto L180;
  }
  L140: m +=  *nblock;
  i__2 =  *nperm;
  i__1 =  *nblock;
  for (i__ = m; i__1 < 0 ? i__ >= i__2: i__ <= i__2; i__ += i__1)
  {
    (*iovect)(n, nblock, &p[p_offset], &i__, &c__1);
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ -  *nblock + j;
      i__4 =  *nperm;
      for (k = 1; k <= i__4; ++k)
      {
        saxpy_(n, &hv[l + k * hv_dim1], &p[j *p_dim1 + 1], &c__1, &vec[k
          *vec_dim1 + 1], &c__1);
        /* L150: */
      }
      /* L160: */
    }
    /* L170: */
  }

  /* ------------------------------------------------------------------ */

  /* THIS SECTION COMPUTES THE RAYLEIGH QUOTIENTS (IN VAL(*,1)) */
  /* AND RESIDUAL NORMS (IN VAL(*,2)) OF THE EIGENVECTORS. */

  L180: m =  *nperm %  *nblock;
  if (m == 0)
  {
    goto L210;
  }
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    scopy_(n, &vec[i__ *vec_dim1 + 1], &c__1, &p[i__ *p_dim1 + 1], &c__1);
    /* L190: */
  }
  (*op)(n, &m, &p[p_offset], &q[q_offset]);
  ++(*nop);
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    val[i__ + val_dim1] = sdot_(n, &p[i__ *p_dim1 + 1], &c__1, &q[i__ *q_dim1 +
      1], &c__1);
    r__1 =  - val[i__ + val_dim1];
    saxpy_(n, &r__1, &p[i__ *p_dim1 + 1], &c__1, &q[i__ *q_dim1 + 1], &c__1);
    val[i__ + (val_dim1 << 1)] = snrm2_(n, &q[i__ *q_dim1 + 1], &c__1);
    /* L200: */
  }
  if (*nperm <  *nblock)
  {
    goto L250;
  }
  L210: ++m;
  i__1 =  *nperm;
  i__2 =  *nblock;
  for (i__ = m; i__2 < 0 ? i__ >= i__1: i__ <= i__1; i__ += i__2)
  {
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ - 1+j;
      scopy_(n, &vec[l *vec_dim1 + 1], &c__1, &p[j *p_dim1 + 1], &c__1);
      /* L220: */
    }
    (*op)(n, nblock, &p[p_offset], &q[q_offset]);
    ++(*nop);
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ - 1+j;
      val[l + val_dim1] = sdot_(n, &p[j *p_dim1 + 1], &c__1, &q[j *q_dim1 + 1],
        &c__1);
      r__1 =  - val[l + val_dim1];
      saxpy_(n, &r__1, &p[j *p_dim1 + 1], &c__1, &q[j *q_dim1 + 1], &c__1);
      val[l + (val_dim1 << 1)] = snrm2_(n, &q[j *q_dim1 + 1], &c__1);
      /* L230: */
    }
    /* L240: */
  }

  /* THIS COMPUTES THE ACCURACY ESTIMATES.  IF AN EIGENVALUE'S RESIDUAL */
  /* NORM INTERVAL LIES INSIDE THE EXCLUDED INTERVAL, THE EIGENPAIR */
  /* IS DELETED.  IF AN EIGENVALUE IS IN THE EXCLUDED INTERVAL BUT ITS */
  /* RESIDUAL NORM INTERVAL OVERLAPS THE BOUNDARY THEN THE EIGENVALUE */
  /* IS SET EQUAL TO THE BOUNDARY, THE RESIDUAL NORM IS RECOMPUTED AND */
  /* IS MADE NEGATIVE, AND 10 IS SUBTRACTED FROM IERR.  A DO LOOP WAS */
  /* NOT USED TO ACCOMODATE THE POSSIBLE DELETION OF EIGENPAIRS. */

  L250: i__ = 0;
  L260: ++i__;
  if (i__ >  *nperm)
  {
    return 0;
  }
  L270: if (val[i__ + val_dim1] <=  *xl || val[i__ + val_dim1] >=  *xr)
  {
    goto L310;
  }
  if (val[i__ + val_dim1] >=  *deltar)
  {
    goto L280;
  }
  if (val[i__ + val_dim1] - val[i__ + (val_dim1 << 1)] >  *xl)
  {
    goto L290;
  }
  temp = val[i__ + val_dim1] -  *xl;
  val[i__ + val_dim1] =  *xl;
  val[i__ + (val_dim1 << 1)] =  - sqrt(temp *temp + val[i__ + (val_dim1 << 1)]
    *val[i__ + (val_dim1 << 1)]);
  *ierr +=  - 10;
  goto L310;

  L280: if (val[i__ + val_dim1] + val[i__ + (val_dim1 << 1)] <  *xr)
  {
    goto L290;
  }
  temp =  *xr - val[i__ + val_dim1];
  val[i__ + val_dim1] =  *xr;
  val[i__ + (val_dim1 << 1)] =  - sqrt(temp *temp + val[i__ + (val_dim1 << 1)]
    *val[i__ + (val_dim1 << 1)]);
  *ierr +=  - 10;
  goto L310;

  L290: --(*nperm);
  /* Computing MIN */
  r__1 =  *deltal, r__2 = val[i__ + val_dim1];
  *deltal = min(r__1, r__2);
  /* Computing MAX */
  r__1 =  *deltar, r__2 = val[i__ + val_dim1];
  *deltar = max(r__1, r__2);
  if (i__ >  *nperm)
  {
    return 0;
  }
  i__2 =  *nperm;
  for (j = i__; j <= i__2; ++j)
  {
    scopy_(n, &vec[(j + 1) *vec_dim1 + 1], &c__1, &vec[j *vec_dim1 + 1],  &c__1)
      ;
    val[j + val_dim1] = val[j + 1+val_dim1];
    val[j + (val_dim1 << 1)] = val[j + 1+(val_dim1 << 1)];
    /* L300: */
  }
  goto L270;

  L310: 
  /* Computing MAX */
  r__1 =  *deltal - val[i__ + val_dim1], r__2 = val[i__ + val_dim1] -  *deltar;
  temp = max(r__1, r__2);
  val[i__ + (val_dim1 << 2)] = (float)0.;
  if (temp > (float)0.)
  {
    val[i__ + (val_dim1 << 2)] = (r__1 = val[i__ + (val_dim1 << 1)], abs(r__1))
      / temp;
  }
  val[i__ + val_dim1 * 3] = val[i__ + (val_dim1 << 2)]*(r__1 = val[i__ + 
    (val_dim1 << 1)], abs(r__1));
  goto L260;

} /* sippla_ */
